clear all
/*Insert directory here*/






/* 1) Create choice sample and run heat estimation */







use ACS1, clear
keep if year==2000 | year==2010 | year==2019
drop if year==2000 & year-builtyr>10
drop if year==2010  & year-builtyr>9
drop if year==2019  & year-builtyr>8
gen furnyear=builtyr
rename year censusyear


/* a) Bring in outside data */


merge m:1 furnyear censusyear statefip geoid using HDDfurngeoid34-2019
keep if _merge==3
drop _merge

merge m:1 furnyear statefip using fuelpricesfurn1934-2019
keep if _merge==3
drop _merge

merge m:1 furnyear statefip censusyear geoid using incfactorfurn2019
keep if _merge==3
drop _merge
gen hhincomefurn=hhincome*incfactorfurn

merge m:1 furnyear using CPIfurn34-2019
keep if _merge==3
drop _merge

local varlist hhincomefurn pnatgasfurn poilfurn pelecfurn ppropanefurn pnatgascounterfurn
foreach a of local varlist{
replace `a'=`a'*255.7/CPI
}

replace hdd=hdd/1000
gen hdd2=hdd*hdd
gen hddfurn2=hddfurn*hddfurn

/* b) create units in structure dummies*/

gen struct2t4=0
replace struct2t4=1  if unitsstr==5 | unitsstr==6
gen struct5=0
replace struct5=1  if unitsstr>=7
gen owner=0
replace owner=1 if ownershp==1
drop ownershp unitsstr

/* e) create household size dummies*/

gen housemem2=0
replace housemem2=1 if numprec==2
gen housemem3=0
replace housemem3=1 if numprec==3
gen housemem4=0
replace housemem4=1 if numprec==4
gen housemem5=0
replace housemem5=1 if numprec==5
gen housemem6=0
replace housemem6=1 if numprec>5

/* f) HDD,rooms,members interaction terms*/

quietly sum hddfurn
gen pnatgascounterfurnxhddfurn=pnatgascounterfurn*(hddfurn-r(mean))
gen poilfurnxhddfurn=poilfurn*(hddfurn-r(mean))
gen pelecfurnxhddfurn=pelecfurn*(hddfurn-r(mean))
gen ppropanefurnxhddfurn=ppropanefurn*(hddfurn-r(mean))

quietly sum rooms
gen pnatgascounterfurnxrooms=pnatgascounterfurn*(rooms-r(mean))
gen poilfurnxrooms=poilfurn*(rooms-r(mean))
gen pelecfurnxrooms=pelecfurn*(rooms-r(mean))
gen ppropanefurnxrooms=ppropanefurn*(rooms-r(mean))

quietly sum numprec
gen pnatgascounterfurnxnumprec=pnatgascounterfurn*(numprec-r(mean))
gen poilfurnxnumprec=poilfurn*(numprec-r(mean))
gen pelecfurnxnumprec=pelecfurn*(numprec-r(mean))
gen ppropanefurnxnumprec=ppropanefurn*(numprec-r(mean))

gen systemchoice=.
replace systemchoice=1 if fuelheat==2
replace systemchoice=2 if fuelheat==5
replace systemchoice=3 if fuelheat==4
replace systemchoice=4 if fuelheat==3
save choiceACS, replace


clear all

program mylogit
  args lnf theta1 theta2 theta3 theta4 
  	 quietly replace `lnf' = ln(exp(`theta1')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==1
     quietly replace `lnf' = ln(exp(`theta2')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==2
	 quietly replace `lnf' = ln(exp(`theta3')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==3
	 	 quietly replace `lnf' = ln(exp(`theta4')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==4
end

use choicecensus, clear

constraint 1 pnatgascounterfurn=poilfurn
constraint 2 pnatgascounterfurnxhddfurn=poilfurnxhddfurn
constraint 3 pnatgascounterfurnxrooms=poilfurnxrooms
constraint 4 pnatgascounterfurnxnumprec=poilfurnxnumprec
constraint 5 ppropanefurn=poilfurn
constraint 6 ppropanefurnxhddfurn=poilfurnxhddfurn
constraint 7 ppropanefurnxrooms=poilfurnxrooms
constraint 8 ppropanefurnxnumprec=poilfurnxnumprec

local covariates hddfurn* housemem* owner rooms hhincomefurn struct*

ml model lf mylogit (gas: systemchoice=pnatgascounterfurn pnatgascounterfurnxhddfurn pnatgascounterfurnxrooms pnatgascounterfurnxnumprec `covariates')(oil: poilfurn poilfurnxhddfurn poilfurnxrooms poilfurnxnumprec, noconstant)(electricity: pelecfurn pelecfurnxhddfurn pelecfurnxrooms pelecfurnxnumprec `covariates')(propane: ppropanefurn ppropanefurnxhddfurn ppropanefurnxrooms ppropanefurnxnumprec `covariates') [pweight=hhwt], constraints (1/8)
ml check
ml search
ml maximize
clear

forvalues y=2001/2013{
use ACS1, clear
rename year censusyear
keep if censusyear==`y'
gen year=`y'
drop if builtyr2==10 & year<2005


/* a) Estimate furnace and water heater replacement years */


gen furnyear=.
gen wateryear=.
gen unif=runiform()

replace furnyear=builtyr if `y'-builtyr<17
replace furnyear=`y' if `y'-builtyr>=17 & unif>=(0/17) & unif<(1/17)
replace furnyear=`y'-1 if `y'-builtyr>=17 & unif>=(1/17) & unif<(2/17)
replace furnyear=`y'-2 if `y'-builtyr>=17 & unif>=(2/17) & unif<(3/17)
replace furnyear=`y'-3 if `y'-builtyr>=17 & unif>=(3/17) & unif<(4/17)
replace furnyear=`y'-4 if `y'-builtyr>=17 & unif>=(4/17) & unif<(5/17)
replace furnyear=`y'-5 if `y'-builtyr>=17 & unif>=(5/17) & unif<(6/17)
replace furnyear=`y'-6 if `y'-builtyr>=17 & unif>=(6/17) & unif<(7/17)
replace furnyear=`y'-7 if `y'-builtyr>=17 & unif>=(7/17) & unif<(8/17)
replace furnyear=`y'-8 if `y'-builtyr>=17 & unif>=(8/17) & unif<(9/17)
replace furnyear=`y'-9 if `y'-builtyr>=17 & unif>=(9/17) & unif<(10/17)
replace furnyear=`y'-10 if `y'-builtyr>=17 & unif>=(10/17) & unif<(11/17)
replace furnyear=`y'-11 if `y'-builtyr>=17 & unif>=(11/17) & unif<(12/17)
replace furnyear=`y'-12 if `y'-builtyr>=17 & unif>=(12/17) & unif<(13/17)
replace furnyear=`y'-13 if `y'-builtyr>=17 & unif>=(13/17) & unif<(14/17)
replace furnyear=`y'-14 if `y'-builtyr>=17 & unif>=(14/17) & unif<(15/17)
replace furnyear=`y'-15 if `y'-builtyr>=17 & unif>=(15/17) & unif<(16/17)
replace furnyear=`y'-16 if `y'-builtyr>=17 & unif>=(16/17) & unif<=(17/17)

replace wateryear=builtyr if `y'-builtyr<10
replace wateryear=`y' if `y'-builtyr>=10 & unif>=(0/10) & unif<(1/10)
replace wateryear=`y'-1 if `y'-builtyr>=10 & unif>=(1/10) & unif<(2/10)
replace wateryear=`y'-2 if `y'-builtyr>=10 & unif>=(2/10) & unif<(3/10)
replace wateryear=`y'-3 if `y'-builtyr>=10 & unif>=(3/10) & unif<(4/10)
replace wateryear=`y'-4 if `y'-builtyr>=10 & unif>=(4/10) & unif<(5/10)
replace wateryear=`y'-5 if `y'-builtyr>=10 & unif>=(5/10) & unif<(6/10)
replace wateryear=`y'-6 if `y'-builtyr>=10 & unif>=(6/10) & unif<(7/10)
replace wateryear=`y'-7 if `y'-builtyr>=10 & unif>=(7/10) & unif<(8/10)
replace wateryear=`y'-8 if `y'-builtyr>=10 & unif>=(8/10) & unif<(9/10)
replace wateryear=`y'-9 if `y'-builtyr>=10 & unif>=(9/10) & unif<=(10/10)
drop unif


/* b) Bring in outside data */


/*current year*/
merge m:1 year statefip geoid using HDDgeoid60-2019
keep if _merge==3
drop _merge

merge m:1 year statefip using fuelprices1960-19
keep if _merge==3
drop _merge

merge m:1 year using CPI34-2019
keep if _merge==3
drop _merge


/*furnace year*/
merge m:1 furnyear censusyear statefip geoid using HDDfurngeoid34-2019
keep if _merge==3
drop _merge

merge m:1 furnyear statefip censusyear geoid using incfactorfurn2019
keep if _merge==3
drop _merge
gen hhincomefurn=hhincome*incfactorfurn

merge m:1 furnyear statefip using fuelpricesfurn1934-2019
keep if _merge==3
drop _merge

merge m:1 furnyear using CPIfurn34-2019
keep if _merge==3
drop _merge


/*water year*/
merge m:1 wateryear censusyear statefip geoid using HDDwatergeoid34-2019
keep if _merge==3
drop _merge

merge m:1 wateryear statefip censusyear geoid using incfactorwater2019
keep if _merge==3
drop _merge
gen hhincomewater=hhincome*incfactorwater

merge m:1 wateryear statefip using fuelpriceswater1934-2019
keep if _merge==3
drop _merge

merge m:1 wateryear using CPIwater34-2019
keep if _merge==3
drop _merge

local varlist hhincome pnatgas poil pelec ppropane pnatgascounter
foreach b of local varlist{
replace `b'=`b'*255.7/CPI
}

local varlist2 hhincomefurn pnatgasfurn poilfurn pelecfurn ppropanefurn pnatgascounterfurn
foreach c of local varlist2{
replace `c'=`c'*255.7/CPIfurn
}

local varlist3 hhincomewater pnatgaswater poilwater pelecwater ppropanewater pnatgascounterwater
foreach d of local varlist3{
replace `d'=`d'*255.7/CPIwater
}

replace hdd=hdd/1000
replace hddfurn=hddfurn/1000
replace hddwater=hddwater/1000
gen hdd2=hdd*hdd
gen hddfurn2=hddfurn*hddfurn
gen hddwater2=hddwater*hddwater
replace cdd=cdd/1000
gen cdd2=cdd*cdd


/* c) create units in structure dummies*/


gen struct2t4=0
replace struct2t4=1  if unitsstr==5 | unitsstr==6
gen struct5=0
replace struct5=1  if unitsstr>=7
gen owner=0
replace owner=1 if ownershp==1
drop ownershp unitsstr


/* d) create age of structure dummies*/


gen ageb1950s=0
gen ageb1960s=0
gen ageb1970s=0
gen ageb1980s=0
gen ageb1990s=0
gen ageb2000s=0
gen ageb2010s=0
replace ageb2010s=1 if builtyr2>=15
replace ageb2000s=1 if builtyr2>=9 & builtyr2<15
replace ageb1990s=1 if builtyr2<=8 & builtyr2>6
replace ageb1980s=1 if builtyr2==6
replace ageb1970s=1 if builtyr2==5
replace ageb1960s=1 if builtyr2==4
replace ageb1950s=1 if builtyr2==3


/* e) create household size dummies*/


gen housemem2=0
replace housemem2=1 if numprec==2
gen housemem3=0
replace housemem3=1 if numprec==3
gen housemem4=0
replace housemem4=1 if numprec==4
gen housemem5=0
replace housemem5=1 if numprec==5
gen housemem6=0
replace housemem6=1 if numprec>5


/* f) HDD,rooms,members interaction terms*/


local varlist furn water
foreach i of local varlist{
quietly sum hdd`i'
gen pnatgascounter`i'xhdd`i'=pnatgascounter`i'*(hdd`i'-r(mean))
gen poil`i'xhdd`i'=poil`i'*(hdd`i'-r(mean))
gen pelec`i'xhdd`i'=pelec`i'*(hdd`i'-r(mean))
gen ppropane`i'xhdd`i'=ppropane`i'*(hdd`i'-r(mean))

quietly sum rooms
gen pnatgascounter`i'xrooms=pnatgascounter`i'*(rooms-r(mean))
gen poil`i'xrooms=poil`i'*(rooms-r(mean))
gen pelec`i'xrooms=pelec`i'*(rooms-r(mean))
gen ppropane`i'xrooms=ppropane`i'*(rooms-r(mean))

quietly sum numprec
gen pnatgascounter`i'xnumprec=pnatgascounter`i'*(numprec-r(mean))
gen poil`i'xnumprec=poil`i'*(numprec-r(mean))
gen pelec`i'xnumprec=pelec`i'*(numprec-r(mean))
gen ppropane`i'xnumprec=ppropane`i'*(numprec-r(mean))
}


/* g) Create predicted values of space heating fuel choice */


gen predg=[gas]_b[pnatgascounterfurn]*pnatgascounterfurn+[gas]_b[pnatgascounterfurnxhddfurn]*pnatgascounterfurnxhddfurn+[gas]_b[pnatgascounterfurnxrooms]*pnatgascounterfurnxrooms+[gas]_b[pnatgascounterfurnxnumprec]*pnatgascounterfurnxnumprec+[gas]_b[hddfurn]*hddfurn+[gas]_b[hddfurn2]*hddfurn2+[gas]_b[housemem2]*housemem2+[gas]_b[housemem3]*housemem3+[gas]_b[housemem4]*housemem4+[gas]_b[housemem5]*housemem5+[gas]_b[housemem6]*housemem6+[gas]_b[owner]*owner+[gas]_b[hhincome]*hhincomefurn+[gas]_b[rooms]*rooms+[gas]_b[struct2t4]*struct2t4+[gas]_b[struct5]*struct5+[gas]_b[_cons]

gen predo=[oil]_b[poilfurn]*poilfurn+[oil]_b[poilfurnxhddfurn]*poilfurnxhddfurn+[oil]_b[poilfurnxrooms]*poilfurnxrooms+[oil]_b[poilfurnxnumprec]*poilfurnxnumprec

gen prede=[electricity]_b[pelecfurn]*pelecfurn+[electricity]_b[pelecfurnxhddfurn]*pelecfurnxhddfurn+[electricity]_b[pelecfurnxrooms]*pelecfurnxrooms+[electricity]_b[pelecfurnxnumprec]*pelecfurnxnumprec+[electricity]_b[hddfurn]*hddfurn+[electricity]_b[hddfurn2]*hddfurn2+[electricity]_b[housemem2]*housemem2+[electricity]_b[housemem3]*housemem3+[electricity]_b[housemem4]*housemem4+[electricity]_b[housemem5]*housemem5+[electricity]_b[housemem6]*housemem6+[electricity]_b[owner]*owner+[electricity]_b[hhincome]*hhincomefurn+[electricity]_b[rooms]*rooms+[electricity]_b[struct2t4]*struct2t4+[electricity]_b[struct5]*struct5+[electricity]_b[_cons]

gen predp=[propane]_b[ppropanefurn]*ppropanefurn+[propane]_b[ppropanefurnxhddfurn]*ppropanefurnxhddfurn+[propane]_b[ppropanefurnxrooms]*ppropanefurnxrooms+[propane]_b[ppropanefurnxnumprec]*ppropanefurnxnumprec+[propane]_b[hddfurn]*hddfurn+[propane]_b[hddfurn2]*hddfurn2+[propane]_b[housemem2]*housemem2+[propane]_b[housemem3]*housemem3+[propane]_b[housemem4]*housemem4+[propane]_b[housemem5]*housemem5+[propane]_b[housemem6]*housemem6+[propane]_b[owner]*owner+[propane]_b[hhincome]*hhincomefurn+[propane]_b[rooms]*rooms+[propane]_b[struct2t4]*struct2t4+[propane]_b[struct5]*struct5+[propane]_b[_cons]

gen choiceg=exp(predg)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choiceo=exp(predo)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choicee=exp(prede)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choicep=exp(predp)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
save dataset`y', replace
}




clear all
/*Insert directory here*/

program mylogit
  args lnf theta1 theta2 theta3 theta4 
  	 quietly replace `lnf' = ln(exp(`theta1')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==1
     quietly replace `lnf' = ln(exp(`theta2')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==2
	 quietly replace `lnf' = ln(exp(`theta3')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==3
	 	 quietly replace `lnf' = ln(exp(`theta4')/(exp(`theta1')+exp(`theta2')+exp(`theta3')+exp(`theta4'))) if $ML_y1==4
end

use choiceACS, clear

constraint 1 pnatgascounterfurn=poilfurn
constraint 2 pnatgascounterfurnxhddfurn=poilfurnxhddfurn
constraint 3 pnatgascounterfurnxrooms=poilfurnxrooms
constraint 4 pnatgascounterfurnxnumprec=poilfurnxnumprec
constraint 5 ppropanefurn=poilfurn
constraint 6 ppropanefurnxhddfurn=poilfurnxhddfurn
constraint 7 ppropanefurnxrooms=poilfurnxrooms
constraint 8 ppropanefurnxnumprec=poilfurnxnumprec

local covariates hddfurn* housemem* owner rooms hhincomefurn struct*

ml model lf mylogit (gas: systemchoice=pnatgascounterfurn pnatgascounterfurnxhddfurn pnatgascounterfurnxrooms pnatgascounterfurnxnumprec `covariates')(oil: poilfurn poilfurnxhddfurn poilfurnxrooms poilfurnxnumprec, noconstant)(electricity: pelecfurn pelecfurnxhddfurn pelecfurnxrooms pelecfurnxnumprec `covariates')(propane: ppropanefurn ppropanefurnxhddfurn ppropanefurnxrooms ppropanefurnxnumprec `covariates') [pweight=hhwt], constraints (1/8)
ml check
ml search
ml maximize
clear

forvalues y=2014/2019{
use ACS1, clear
rename year censusyear
keep if censusyear==`y'
gen year=`y'
drop if builtyr2==10 & year<2005


/* a) Estimate furnace and water heater replacement years */


gen furnyear=.
gen wateryear=.
gen unif=runiform()

replace furnyear=builtyr if `y'-builtyr<17
replace furnyear=`y' if `y'-builtyr>=17 & unif>=(0/17) & unif<(1/17)
replace furnyear=`y'-1 if `y'-builtyr>=17 & unif>=(1/17) & unif<(2/17)
replace furnyear=`y'-2 if `y'-builtyr>=17 & unif>=(2/17) & unif<(3/17)
replace furnyear=`y'-3 if `y'-builtyr>=17 & unif>=(3/17) & unif<(4/17)
replace furnyear=`y'-4 if `y'-builtyr>=17 & unif>=(4/17) & unif<(5/17)
replace furnyear=`y'-5 if `y'-builtyr>=17 & unif>=(5/17) & unif<(6/17)
replace furnyear=`y'-6 if `y'-builtyr>=17 & unif>=(6/17) & unif<(7/17)
replace furnyear=`y'-7 if `y'-builtyr>=17 & unif>=(7/17) & unif<(8/17)
replace furnyear=`y'-8 if `y'-builtyr>=17 & unif>=(8/17) & unif<(9/17)
replace furnyear=`y'-9 if `y'-builtyr>=17 & unif>=(9/17) & unif<(10/17)
replace furnyear=`y'-10 if `y'-builtyr>=17 & unif>=(10/17) & unif<(11/17)
replace furnyear=`y'-11 if `y'-builtyr>=17 & unif>=(11/17) & unif<(12/17)
replace furnyear=`y'-12 if `y'-builtyr>=17 & unif>=(12/17) & unif<(13/17)
replace furnyear=`y'-13 if `y'-builtyr>=17 & unif>=(13/17) & unif<(14/17)
replace furnyear=`y'-14 if `y'-builtyr>=17 & unif>=(14/17) & unif<(15/17)
replace furnyear=`y'-15 if `y'-builtyr>=17 & unif>=(15/17) & unif<(16/17)
replace furnyear=`y'-16 if `y'-builtyr>=17 & unif>=(16/17) & unif<=(17/17)

replace wateryear=builtyr if `y'-builtyr<10
replace wateryear=`y' if `y'-builtyr>=10 & unif>=(0/10) & unif<(1/10)
replace wateryear=`y'-1 if `y'-builtyr>=10 & unif>=(1/10) & unif<(2/10)
replace wateryear=`y'-2 if `y'-builtyr>=10 & unif>=(2/10) & unif<(3/10)
replace wateryear=`y'-3 if `y'-builtyr>=10 & unif>=(3/10) & unif<(4/10)
replace wateryear=`y'-4 if `y'-builtyr>=10 & unif>=(4/10) & unif<(5/10)
replace wateryear=`y'-5 if `y'-builtyr>=10 & unif>=(5/10) & unif<(6/10)
replace wateryear=`y'-6 if `y'-builtyr>=10 & unif>=(6/10) & unif<(7/10)
replace wateryear=`y'-7 if `y'-builtyr>=10 & unif>=(7/10) & unif<(8/10)
replace wateryear=`y'-8 if `y'-builtyr>=10 & unif>=(8/10) & unif<(9/10)
replace wateryear=`y'-9 if `y'-builtyr>=10 & unif>=(9/10) & unif<=(10/10)
drop unif


/* b) Bring in outside data */


/*current year*/
merge m:1 year statefip geoid using HDDgeoid60-2019
keep if _merge==3
drop _merge

merge m:1 year statefip using fuelprices1960-19
keep if _merge==3
drop _merge

merge m:1 year using CPI34-2019
keep if _merge==3
drop _merge


/*furnace year*/
merge m:1 furnyear censusyear statefip geoid using HDDfurngeoid34-2019
keep if _merge==3
drop _merge

merge m:1 furnyear statefip censusyear geoid using incfactorfurn2019
keep if _merge==3
drop _merge
gen hhincomefurn=hhincome*incfactorfurn

merge m:1 furnyear statefip using fuelpricesfurn1934-2019
keep if _merge==3
drop _merge

merge m:1 furnyear using CPIfurn34-2019
keep if _merge==3
drop _merge


/*water year*/
merge m:1 wateryear censusyear statefip geoid using HDDwatergeoid34-2019
keep if _merge==3
drop _merge

merge m:1 wateryear statefip censusyear geoid using incfactorwater2019
keep if _merge==3
drop _merge
gen hhincomewater=hhincome*incfactorwater

merge m:1 wateryear statefip using fuelpriceswater1934-2019
keep if _merge==3
drop _merge

merge m:1 wateryear using CPIwater34-2019
keep if _merge==3
drop _merge

local varlist hhincome pnatgas poil pelec ppropane pnatgascounter
foreach b of local varlist{
replace `b'=`b'*255.7/CPI
}

local varlist2 hhincomefurn pnatgasfurn poilfurn pelecfurn ppropanefurn pnatgascounterfurn
foreach c of local varlist2{
replace `c'=`c'*255.7/CPIfurn
}

local varlist3 hhincomewater pnatgaswater poilwater pelecwater ppropanewater pnatgascounterwater
foreach d of local varlist3{
replace `d'=`d'*255.7/CPIwater
}

replace hdd=hdd/1000
replace hddfurn=hddfurn/1000
replace hddwater=hddwater/1000
gen hdd2=hdd*hdd
gen hddfurn2=hddfurn*hddfurn
gen hddwater2=hddwater*hddwater
replace cdd=cdd/1000
gen cdd2=cdd*cdd


/* c) create units in structure dummies*/


gen struct2t4=0
replace struct2t4=1  if unitsstr==5 | unitsstr==6
gen struct5=0
replace struct5=1  if unitsstr>=7
gen owner=0
replace owner=1 if ownershp==1
drop ownershp unitsstr


/* d) create age of structure dummies*/


gen ageb1950s=0
gen ageb1960s=0
gen ageb1970s=0
gen ageb1980s=0
gen ageb1990s=0
gen ageb2000s=0
gen ageb2010s=0
replace ageb2010s=1 if builtyr2>=15
replace ageb2000s=1 if builtyr2>=9 & builtyr2<15
replace ageb1990s=1 if builtyr2<=8 & builtyr2>6
replace ageb1980s=1 if builtyr2==6
replace ageb1970s=1 if builtyr2==5
replace ageb1960s=1 if builtyr2==4
replace ageb1950s=1 if builtyr2==3


/* e) create household size dummies*/


gen housemem2=0
replace housemem2=1 if numprec==2
gen housemem3=0
replace housemem3=1 if numprec==3
gen housemem4=0
replace housemem4=1 if numprec==4
gen housemem5=0
replace housemem5=1 if numprec==5
gen housemem6=0
replace housemem6=1 if numprec>5


/* f) HDD,rooms,members interaction terms*/


local varlist furn water
foreach i of local varlist{
quietly sum hdd`i'
gen pnatgascounter`i'xhdd`i'=pnatgascounter`i'*(hdd`i'-r(mean))
gen poil`i'xhdd`i'=poil`i'*(hdd`i'-r(mean))
gen pelec`i'xhdd`i'=pelec`i'*(hdd`i'-r(mean))
gen ppropane`i'xhdd`i'=ppropane`i'*(hdd`i'-r(mean))

quietly sum rooms
gen pnatgascounter`i'xrooms=pnatgascounter`i'*(rooms-r(mean))
gen poil`i'xrooms=poil`i'*(rooms-r(mean))
gen pelec`i'xrooms=pelec`i'*(rooms-r(mean))
gen ppropane`i'xrooms=ppropane`i'*(rooms-r(mean))

quietly sum numprec
gen pnatgascounter`i'xnumprec=pnatgascounter`i'*(numprec-r(mean))
gen poil`i'xnumprec=poil`i'*(numprec-r(mean))
gen pelec`i'xnumprec=pelec`i'*(numprec-r(mean))
gen ppropane`i'xnumprec=ppropane`i'*(numprec-r(mean))
}


/* g) Create predicted values of space heating fuel choice */


gen predg=[gas]_b[pnatgascounterfurn]*pnatgascounterfurn+[gas]_b[pnatgascounterfurnxhddfurn]*pnatgascounterfurnxhddfurn+[gas]_b[pnatgascounterfurnxrooms]*pnatgascounterfurnxrooms+[gas]_b[pnatgascounterfurnxnumprec]*pnatgascounterfurnxnumprec+[gas]_b[hddfurn]*hddfurn+[gas]_b[hddfurn2]*hddfurn2+[gas]_b[housemem2]*housemem2+[gas]_b[housemem3]*housemem3+[gas]_b[housemem4]*housemem4+[gas]_b[housemem5]*housemem5+[gas]_b[housemem6]*housemem6+[gas]_b[owner]*owner+[gas]_b[hhincome]*hhincomefurn+[gas]_b[rooms]*rooms+[gas]_b[struct2t4]*struct2t4+[gas]_b[struct5]*struct5+[gas]_b[_cons]

gen predo=[oil]_b[poilfurn]*poilfurn+[oil]_b[poilfurnxhddfurn]*poilfurnxhddfurn+[oil]_b[poilfurnxrooms]*poilfurnxrooms+[oil]_b[poilfurnxnumprec]*poilfurnxnumprec

gen prede=[electricity]_b[pelecfurn]*pelecfurn+[electricity]_b[pelecfurnxhddfurn]*pelecfurnxhddfurn+[electricity]_b[pelecfurnxrooms]*pelecfurnxrooms+[electricity]_b[pelecfurnxnumprec]*pelecfurnxnumprec+[electricity]_b[hddfurn]*hddfurn+[electricity]_b[hddfurn2]*hddfurn2+[electricity]_b[housemem2]*housemem2+[electricity]_b[housemem3]*housemem3+[electricity]_b[housemem4]*housemem4+[electricity]_b[housemem5]*housemem5+[electricity]_b[housemem6]*housemem6+[electricity]_b[owner]*owner+[electricity]_b[hhincome]*hhincomefurn+[electricity]_b[rooms]*rooms+[electricity]_b[struct2t4]*struct2t4+[electricity]_b[struct5]*struct5+[electricity]_b[_cons]

gen predp=[propane]_b[ppropanefurn]*ppropanefurn+[propane]_b[ppropanefurnxhddfurn]*ppropanefurnxhddfurn+[propane]_b[ppropanefurnxrooms]*ppropanefurnxrooms+[propane]_b[ppropanefurnxnumprec]*ppropanefurnxnumprec+[propane]_b[hddfurn]*hddfurn+[propane]_b[hddfurn2]*hddfurn2+[propane]_b[housemem2]*housemem2+[propane]_b[housemem3]*housemem3+[propane]_b[housemem4]*housemem4+[propane]_b[housemem5]*housemem5+[propane]_b[housemem6]*housemem6+[propane]_b[owner]*owner+[propane]_b[hhincome]*hhincomefurn+[propane]_b[rooms]*rooms+[propane]_b[struct2t4]*struct2t4+[propane]_b[struct5]*struct5+[propane]_b[_cons]

gen choiceg=exp(predg)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choiceo=exp(predo)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choicee=exp(prede)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
gen choicep=exp(predp)/(exp(predg)+exp(predo)+exp(prede)+exp(predp))
save dataset`y', replace
}






/* 2) Run RECS water heater choice model for ACS */






clear all
/*Insert directory here*/

program mylogit
  args lnf theta1 theta2
  	 quietly replace `lnf' = ln(exp(`theta1')/(exp(`theta1')+exp(`theta2'))) if $ML_y1==1
     quietly replace `lnf' = ln(exp(`theta2')/(exp(`theta1')+exp(`theta2'))) if $ML_y1==2
end


/* a) Water heater with gas furnace */


use RECSadj, clear
keep if fuelheat==1
keep if fuelh2o==1 | fuelh2o==5

gen systemchoice=.
replace systemchoice=1 if fuelh2o==5
replace systemchoice=2 if fuelh2o==1

local covariates hdd* owner rooms hhincome struct* housemem*

ml model lf mylogit (elec: systemchoice=elecprice elecpricexhdd elecpricexrooms elecpricexhhmem  `covariates')(gas: gasprice gaspricexhdd gaspricexrooms gaspricexhhmem, noconstant) [pweight=nweight]
ml check
ml search
ml maximize
clear

forvalues y=2001/2019{
use dataset`y', clear

/* electric with gas furnace */

gen predwaterge=[elec]_b[elecprice]*pelecwater+[elec]_b[elecpricexhdd]*pelecwaterxhddwater+[elec]_b[elecpricexrooms]*pelecwaterxrooms+[elec]_b[elecpricexhhmem]*pelecwaterxnumprec+[elec]_b[hdd]*hddwater+[elec]_b[hdd2]*hddwater2+[elec]_b[housemem2]*housemem2+[elec]_b[housemem3]*housemem3+[elec]_b[housemem4]*housemem4+[elec]_b[housemem5]*housemem5+[elec]_b[housemem6]*housemem6+[elec]_b[owner]*owner+[elec]_b[hhincome]*hhincome+[elec]_b[rooms]*rooms+[elec]_b[struct2t4]*struct2t4+[elec]_b[struct5]*struct5+[elec]_b[_cons]

/* gas with gas furnace */
gen predwatergg=[gas]_b[gasprice]*pnatgascounterwater+[gas]_b[gaspricexhdd]*pnatgascounterwaterxhddwater+[gas]_b[gaspricexrooms]*pnatgascounterwaterxrooms+[gas]_b[gaspricexhhmem]*pnatgascounterwaterxnumprec

gen choicewaterge=exp(predwaterge)/(exp(predwaterge)+exp(predwatergg))
gen choicewatergg=exp(predwatergg)/(exp(predwaterge)+exp(predwatergg))
drop predwaterge predwatergg
save dataset`y', replace
}


/* b) Water heater with electric furnace */

use RECSadj, clear
keep if fuelheat==5
keep if fuelh2o==1 | fuelh2o==5

gen systemchoice=.
replace systemchoice=1 if fuelh2o==1
replace systemchoice=2 if fuelh2o==5

local covariates hdd* owner rooms hhincome struct* housemem*

ml model lf mylogit (gas: systemchoice=gasprice gaspricexhdd gaspricexrooms gaspricexhhmem  `covariates')(elec: elecprice elecpricexhdd elecpricexrooms elecpricexhhmem, noconstant) [pweight=nweight]
ml check
ml search
ml maximize
clear

forvalues y=2001/2019{
use dataset`y', clear

/* gas with electric furnace */
gen predwatereg=[gas]_b[gasprice]*pnatgascounterwater+[gas]_b[gaspricexhdd]*pnatgascounterwaterxhddwater+[gas]_b[gaspricexrooms]*pnatgascounterwaterxrooms+[gas]_b[gaspricexhhmem]*pnatgascounterwaterxnumprec+[gas]_b[hdd]*hddwater+[gas]_b[hdd2]*hddwater2+[gas]_b[housemem2]*housemem2+[gas]_b[housemem3]*housemem3+[gas]_b[housemem4]*housemem4+[gas]_b[housemem5]*housemem5+[gas]_b[housemem6]*housemem6+[gas]_b[owner]*owner+[gas]_b[hhincome]*hhincome+[gas]_b[rooms]*rooms+[gas]_b[struct2t4]*struct2t4+[gas]_b[struct5]*struct5+[gas]_b[_cons]

/*electric with electric furnace*/
gen predwateree=[elec]_b[elecprice]*pelecwater+[elec]_b[elecpricexhdd]*pelecwaterxhddwater+[elec]_b[elecpricexrooms]*pelecwaterxrooms+[elec]_b[elecpricexhhmem]*pelecwaterxnumprec

gen choicewatereg=exp(predwatereg)/(exp(predwatereg)+exp(predwateree))
gen choicewateree=exp(predwateree)/(exp(predwatereg)+exp(predwateree))
drop predwateree predwatereg
save dataset`y', replace
}


/* c) Water heater with oil furnace*/

use RECSadj, clear
keep if fuelheat==3
keep if fuelh2o==3 | fuelh2o==5

gen systemchoice=.
replace systemchoice=1 if fuelh2o==5
replace systemchoice=2 if fuelh2o==3

local covariates hdd* owner rooms hhincome struct* housemem*

ml model lf mylogit (elec: systemchoice=elecprice elecpricexhdd elecpricexrooms elecpricexhhmem  `covariates')(oil: oilprice oilpricexhdd oilpricexrooms oilpricexhhmem, noconstant) [pweight=nweight]
ml check
ml search
ml maximize
clear

forvalues y=2001/2019{
use dataset`y', clear

/*electric with oil furnace*/
gen predwateroe=[elec]_b[elecprice]*pelecwater+[elec]_b[elecpricexhdd]*pelecwaterxhddwater+[elec]_b[elecpricexrooms]*pelecwaterxrooms+[elec]_b[elecpricexhhmem]*pelecwaterxnumprec+[elec]_b[hdd]*hddwater+[elec]_b[hdd2]*hddwater2+[elec]_b[housemem2]*housemem2+[elec]_b[housemem3]*housemem3+[elec]_b[housemem4]*housemem4+[elec]_b[housemem5]*housemem5+[elec]_b[housemem6]*housemem6+[elec]_b[owner]*owner+[elec]_b[hhincome]*hhincome+[elec]_b[rooms]*rooms+[elec]_b[struct2t4]*struct2t4+[elec]_b[struct5]*struct5+[elec]_b[_cons]

/*oil with oil furnace*/
gen predwateroo=[oil]_b[oilprice]*poilwater+[oil]_b[oilpricexhdd]*poilwaterxhddwater+[oil]_b[oilpricexrooms]*poilwaterxrooms+[oil]_b[oilpricexhhmem]*poilwaterxnumprec

gen choicewateroe=exp(predwateroe)/(exp(predwateroe)+exp(predwateroo))
gen choicewateroo=exp(predwateroo)/(exp(predwateroe)+exp(predwateroo))
drop predwateroo predwateroe
save dataset`y', replace
}


/* d) Propane furnace water heater fuel */

use RECSadj, clear
keep if fuelheat==2
keep if fuelh2o==2 | fuelh2o==5

gen systemchoice=.
replace systemchoice=1 if fuelh2o==5
replace systemchoice=2 if fuelh2o==2

local covariates hdd* owner rooms hhincome struct* housemem*

ml model lf mylogit (elec: systemchoice=elecprice elecpricexhdd elecpricexrooms elecpricexhhmem  `covariates')(prop: propprice proppricexhdd proppricexrooms proppricexhhmem, noconstant) [pweight=nweight]
ml check
ml search
ml maximize
clear

forvalues y=2001/2019{
use dataset`y', clear

/*electric with prop furnace*/
gen predwaterpe=[elec]_b[elecprice]*pelecwater+[elec]_b[elecpricexhdd]*pelecwaterxhddwater+[elec]_b[elecpricexrooms]*pelecwaterxrooms+[elec]_b[elecpricexhhmem]*pelecwaterxnumprec+[elec]_b[hdd]*hddwater+[elec]_b[hdd2]*hddwater2+[elec]_b[housemem2]*housemem2+[elec]_b[housemem3]*housemem3+[elec]_b[housemem4]*housemem4+[elec]_b[housemem5]*housemem5+[elec]_b[housemem6]*housemem6+[elec]_b[owner]*owner+[elec]_b[hhincome]*hhincome+[elec]_b[rooms]*rooms+[elec]_b[struct2t4]*struct2t4+[elec]_b[struct5]*struct5+[elec]_b[_cons]

/*prop with prop furnace*/
gen predwaterpp=[prop]_b[propprice]*ppropanewater+[prop]_b[proppricexhdd]*ppropanewaterxhddwater+[prop]_b[proppricexrooms]*ppropanewaterxrooms+[prop]_b[proppricexhhmem]*ppropanewaterxnumprec

gen choicewaterpe=exp(predwaterpe)/(exp(predwaterpe)+exp(predwaterpp))
gen choicewaterpp=exp(predwaterpp)/(exp(predwaterpe)+exp(predwaterpp))
drop predwaterpp predwaterpe

gen choicewaterg=choiceg*choicewatergg+choicee*choicewatereg
gen choicewatero=choiceo*choicewateroo
gen choicewaterp=choicep*choicewaterpp
gen choicewatere=choicee*choicewateree+choiceg*choicewaterge+choiceo*choicewateroe+choicep*choicewaterpe
save dataset`y', replace
}





/* 3) Generate Demand ACS*/






/* a) Electric Space Heating Demand*/

forvalues y=2001/2019{
use RECSadj, clear
gen id1=0
replace id1=1 if year>=1997 & year<=2005
gen id2=0
replace id2=1 if year>=2005 & year<=2009
gen id3=0
replace id3=1 if year>=2009

drop if id1!=1 & `y'>=2001 & `y'<=2005
drop if id2!=1 & `y'>=2006 & `y'<=2011
drop if id3!=1 & `y'>2011 & `y'<=2019
drop id1 id2 id3

keep if fuelheat==5
reg heatcons heatprice hdd hdd2 hhincome rooms build* struct* owner [pweight=nweight], robust

use dataset`y', clear

gen predelecspacecons=_b[heatprice]*pelec+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950s+_b[build1960s]*ageb1960s+_b[build1970s]*ageb1970s+_b[build1980s]*ageb1980s+_b[build1990s]*ageb1990s+_b[build2000s]*ageb2000+_b[_cons]
replace predelecspacecons=predelecspacecons+_b[build2010s]*ageb2010 if `y'>=2010
replace predelecspacecons=0 if predelecspacecons<0

save dataset`y', replace
}


/* b) NatGas & Oil Space Heating Demand*/

forvalues y=2001/2019{
use RECSadj, clear
gen id1=0
replace id1=1 if year>=1997 & year<=2005
gen id2=0
replace id2=1 if year>=1997 & year<=2005
gen id3=0
replace id3=1 if year>=2005
gen id4=0
replace id4=1 if year>=2009

drop if id1!=1 & `y'==2001
drop if id2!=1 & `y'>2001 & `y'<=2008
drop if id3!=1 & `y'>2008 & `y'<=2011
drop if id4!=1 & `y'>2011 & `y'<=2019
drop id1 id2 id3 id4
keep if fuelheat==1 | fuelheat==3
reg heatcons heatprice hdd hdd2 hhincome rooms build* struct* owner [pweight=nweight], robust

use dataset`y', clear

local varlist natgas oil
foreach j of local varlist{
gen pred`j'spacecons=_b[heatprice]*p`j'+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950s+_b[build1960s]*ageb1960s+_b[build1970s]*ageb1970s+_b[build1980s]*ageb1980s+_b[build1990s]*ageb1990s+_b[build2000s]*ageb2000s+_b[_cons]
replace pred`j'spacecons=pred`j'spacecons+_b[build2010s]*ageb2010 if `y'>=2010
replace pred`j'spacecons=0 if pred`j'spacecons<0
}

save dataset`y', replace
}


/* c) Propane Space Heating Demand*/

forvalues y=2001/2019{
use RECSadj, clear
gen id1=0
replace id1=1 if year>=1997 & year<=2005
gen id2=0
replace id2=1 if year>=2001 & year<=2009
gen id3=0
replace id3=1 if year>=2009

drop if id1!=1 & `y'>=2001 & `y'<=2005
drop if id2!=1 & `y'>=2004 & `y'<=2013
drop if id3!=1 & `y'>2013 & `y'<=2019
drop id1 id2 id3
keep if fuelheat==2
reg heatcons heatprice hdd hdd2 hhincome rooms build* struct* owner [pweight=nweight], robust

use dataset`y', clear

gen predpropanespacecons=_b[heatprice]*ppropane+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950s+_b[build1960s]*ageb1960s+_b[build1970s]*ageb1970s+_b[build1980s]*ageb1980s+_b[build1990s]*ageb1990s+_b[build2000s]*ageb2000s+_b[_cons]
replace predpropanespacecons=predpropanespacecons+_b[build2010s]*ageb2010s if `y'>=2010
replace predpropanespacecons=0 if predpropanespacecons<0
save dataset`y', replace
}


/* d) Generate Water Heating Demand */


/*Electric Water Heating*/
forvalues y=2001/2019{
use RECSadj, clear
gen id1=0
replace id1=1 if year>=1997 & year<=2005
gen id2=0
replace id2=1 if year>=2005 & year<=2009
gen id3=0
replace id3=1 if year>=2009

drop if id1!=1 & `y'>=2001 & `y'<=2005
drop if id2!=1 & `y'>=2006 & `y'<=2011
drop if id3!=1 & `y'>2011 & `y'<=2019
drop id1 id2 id3
keep if fuelh2o==5
drop if heatprice<1 | heatprice>300
drop if watercons<0.5 | watercons>80

reg watercons elecprice hdd hdd2 hhincome rooms struct* housemem* owner [pweight=nweight], robust

use dataset`y', clear

gen predelecwatercons=_b[elecprice]*pelec+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[_cons]
replace predelecwatercons=0 if predelecwatercons<0

save dataset`y', replace


/*Fossil Fuel Water Heating*/

use RECSadj, clear
gen id1=0
replace id1=1 if year>=1997 & year<=2005
gen id2=0
replace id2=1 if year>=2005 & year<=2009
gen id3=0
replace id3=1 if year>=2009

drop if id1!=1 & `y'>=2001 & `y'<=2005
drop if id2!=1 & `y'>=2006 & `y'<=2011
drop if id3!=1 & `y'>2011 & `y'<=2019
drop id1 id2 id3
drop if fuelh2o==5
drop if heatprice<1 | heatprice>300
drop if watercons<0.5 | watercons>80

reg watercons waterheatprice hdd hdd2 hhincome rooms struct* housemem* owner [pweight=nweight], robust

use dataset`y', clear

local varlist2 propane oil
foreach j of local varlist2{
gen pred`j'watercons=_b[waterheatprice]*p`j'+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[_cons]
replace pred`j'watercons=0 if pred`j'watercons<0
}

gen prednatgaswatercons=_b[waterheatprice]*pnatgascounter+_b[hdd]*hdd+_b[hdd2]*hdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[_cons]
replace prednatgaswatercons=0 if prednatgaswatercons<0

save dataset`y', replace



/* e) Generate Non-heating electricity demand */
use RECSadj, clear

gen id1=0
replace id1=1 if year<=2001
gen id2=0
replace id2=1 if year>=1993 & year<=2005
gen id3=0
replace id3=1 if year>=2001
gen id4=0
replace id4=1 if year>=2005

drop if id1!=1 & `y'==2001
drop if id2!=1 & `y'>=2002 & `y'<=2004
drop if id3!=1 & `y'>2004 & `y'<=2009
drop if id4!=1 & `y'>2009
drop id1 id2 id3 id4
reg comeleccons elecprice cdd cdd2 hhincome rooms build* struct* owner [pweight=nweight], robust

use dataset`y', clear

gen predcomeleccons=_b[elecprice]*pelec+_b[cdd]*cdd+_b[cdd2]*cdd2+_b[owner]*owner+_b[hhincome]*hhincome+_b[rooms]*rooms+_b[struct2t4]*struct2t4+_b[struct5]*struct5+_b[build1950s]*ageb1950+_b[build1960s]*ageb1960+_b[build1970s]*ageb1970+_b[build1980s]*ageb1980+_b[build1990s]*ageb1990+_b[build2000s]*ageb2000+_b[_cons]
replace predcomeleccons=predcomeleccons+_b[build2010s]*ageb2010 if `y'>=2010
replace predcomeleccons=0 if predcomeleccons<0

save dataset`y', replace
}

